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ABSTRACT 

Direct measurements of the extragalactic background light (EBL) in the near-IR 
to mid-IR waveband are extremely difficult due to an overwhelming foreground from 
the zodiacal light that outshines the faint cosmological diffuse radiation field by more 
than an order of magnitude. Indirect constraints on the EBL are provided by 7-ray 
observations of AGN. Using the combination of the Fermi Gamma-Ray Space Tele- 
scope together with the current generation of ground-based air Cherenkov telescopes 
(H.E.S.S., MAGIC, and VERITAS) provides unprecedented sensitivity and spectral 
coverage for constraining the EBL in the near- to mid-IR. In this paper we present 
new limits on the EBL based on the analysis of the broad-band spectra of a select set 
of 7-ray blazars covering 200 MeV to several TeV. The EBL intensity at 15 //m is con- 
strained to be 1.36 ± 0.58nWm~^ sr~^. We find that the fast evolution and baseline 
EBL models of Stecker et al. (2006), as well as the model of Kneiske et al. (2004), 
predict significantly higher EBL intensities in the mid-IR (15 ^m) than is allowed by 
the constraints derived here. In addition, the model of Franceschini et al. (2008) and 
the fiducial model of Dommguez et al. (2011) predict near- to mid-IR ratios smaller 
than that predicted by our analysis. Namely, their intensities in the near-IR are too low 
while their intensities in the mid-IR arc marginally too high. All of the aforementioned 
models are inconsistent with our analysis at the >3(7 level. 

Subject headings: diffuse radiation — galaxies: active — gamma rays: galaxies — 
infrared: diffuse background 

1. Introduction 

The extragalactic background light (EBL) contains all radiative energy releases, from nuclear 
and accretion processes, that have occurred since the epoch of recombination. While the EBL 
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is the second most dominant cosmological radiation field in our universe, surpassed only by the 
cosmic microwave background (CMB), it has escaped detection through direct measurements at 
most wavelengths. Its spectrum is bimodal with one component peaking at ~1 /im, originating 
from radiation released through the formation of heavy elements and the accretion of matter onto 
black holes in active galactic nuclei (AGNs), and a second component peaking at ~100//m, created 
through the absorption of UV/optical light that is re-radiated by dust at infrared wavelengths. 
The energy spectrum contains, therefore, important information about the cosmic evolution of 
these energy sources, the obscuring dust, and the relative contributions of starburst galaxies and 
AGNs to the energy released over cosmic time. In addition, the EBL contains information on the 
rate of core collapse supernovae which provides an important constraint on the normalization of the 
diffuse supernova neutrino background ( Horiuchi et al.|["2009 >. The trough in the EBL spectrum at 
~15 /im is caused by the decrease in the stellar and AGN UV-optical emission towards the mid-IR, 
and the paucity of hot dust capable of emitting at these wavelengths. This regime is also the most 
difficult for direct measurements of the EBL due to an overwhelming foreground from zodiacal 
light. For detailed reviews of the origin, measurements, modeling, and cosmological implications of 



the EBL see Hauser & Dwek (2001) and Kashlinsky (2005). 



Very high energy (VHE) gamma-rays (~0.1— 10 TeV) interact, via pair production, with photons 
in the near- to mid-IR regime of the EBL (i.e., 7vhe7ebl ~^ e+e") dGould Schreder||l967| >. Due 
to the energy dependence of the 7-7 optical depth, the observed VHE emission spectra of blazars 
are softer than the intrinsic spectra. Information regarding the EBL spectral energy distribution 
(SED) can be gleaned from this spectral absorption using a variety of approaches. A common 
technique employed assumes a theoretically based limit on the hardness of the intrinsic spectrum, 
and hence places a limit on the maximum level of EBL absorption given the softness of the observed 



spectrum (Aharonian et al. 2006, 2007d). This highest allowed level of absorption translates into 



an upper limit on the EBL intensity. 



Another technique, being an expanded approach to that of Aharonian et al. (1999), utilized 



the TeV spectra of nearby blazars Markarian 421 and Markarian 501 {z ~ 0.03) (Dwek & Krennrich 



2005 ) . The authors identified EBL scenarios producing an exponential rise in luminosity with energy 



in absorption-corrected spectra. Such an exponential rise is highly unlikely, based on standard 



synchrotron self-Compton (Maraschi et al. 1992 Bloom Sz Marscher 1996) and external inverse- 



Compton ( Dermer fc Schlickeiser|1993 Sikora et al.|1994 ) models of blazars, and hence the authors 
conclude these scenarios are unrealistic. 

In this paper we use two different approaches to constrain the EBL, one of which only recently 
became feasible as a result of Fermi measurements of blazar spectra in the 0.1 to 100 GeV energy 
regime. At these energies, blazar spectra are typically well characterized by a power-law and 
are a good representations of the intrinsic spectra, since the EBL causes a negligible amount of 
attenuation at these energies. The first approach, referred to hereafter as the spectral shape method, 
or Method 1, assumes that the same power-law characterizing the blazar spectrum at GeV energies 
extends up to TeV energies, while also allowing for curvature that softens the spectrum at very 
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high energies. Consequently, only EBL SEDs yielding intrinsic blazar spectra consistent with the 
TeV-extended Fermi power-law are considered viable. The use of the power-law assumption is, of 
course, only valid over a limited energy regime since the inverse-Compton (IC) component shows 
a turnover close to its peak energy. However, spectra with IC peaks in the multi-TeV regime are 



generally well characterized by a power-law below the peak (Samuelson et al. 1998; Aharonian 



et al. 1999 Krennrich et al. 2001 Aharonian et al. 2002). Similar approaches have been used by 



Georganopoulos et al. (2010) and Mankuzhiyil et al. (2010) to constrain the EBL absorption optical 



depth for VHE photons and by Prandini et al. (2010) to constrain the distances of blazars with 
unknown redshifts. 

The second approach, referred to as the TeV spectral break method, or Method 2, makes a less 
stringent assumption, namely that the intrinsic blazar spectrum at TeV energies is characterized by 
a single power-law, with no restrictions placed on the spectral index. We then search for evidence 
of a break at ~1 TeV in observed blazar spectra, caused by EBL absorption, and compare this with 
the spectral breaks predicted by different EBL scenarios. The EBL SEDs producing blazar spectral 
breaks consistent with observations are considered viable. 

We characterize the family of EBL realizations by the 1.6 /im and 15;um intensities, and the 
ratio between the two. As we show later on, the two approaches are complementary. The spectral 
shape method is most sensitive to the overall EBL intensity, whereas the TeV spectral break 
method is mostly sensitive to the relative difference between the 1.6 ^m and 15 /im intensities, that 
is, the slope of the EBL between these two wavelengths. The two methods therefore probe different 
features of the EBL SED, ruling out different regions of parameter space, providing new constraints 
on the EBL. In Section [2j we discuss the range of EBL parameter space tested. Section [3] describes 
the two analysis methods used and how they can be combined to improve constraints on the EBL. 
The results of our analysis are presented in Section |4j with associated interesting and important 
caveats outlined in Section [5] Section [6] provides some suggestions for improving constraints on the 
EBL with future blazar measurements. Finally, a discussion of our results, along with concluding 
remarks, is provided in Section [7j 



2. EBL Models 



A variety of models exist describing the SED of the EBL (e.g., Kneiske et al. (2002); Primack 
et al. (2005); Stecker et al. (2006); Franceschini et al. ( 2008[ )). The techniques used to generate 
these models fall into three basic categories: forward evolution, backward evolution, and observed 
evolution over a particular redshift range. The forward evolution method (e.g., Primack et al 



(2005)) makes use of early structure formation scenarios to predict the evolution of galactic lumi- 



nosity functions forward in time. Backward evolution techniques (e.g., Stecker et al. (2006)) begin 
with the existing galaxy population and model their luminosity functions backward in time. The 



third approach to modeling uses deep galaxy surveys (e.g., Kneiske et al. (2002); Franceschini et al 



(2008)) or tracers of chemical evolution to infer the cosmic star formation rate and compute the 



- 4 - 



100 



a 




1000 



10 



K 1 



0.1 



0.1 




1 

E [TeV] 



10 



Fig. 1. — Left: EBL intensity versus photon wavelength. The shaded region indicates the range of 
scenarios tested. The thick sohd hne indicates the basehne shape used, from which all other scaled 
shapes are generated. For clarity, two additional models are shown (dotted and dashed) illustrating 
the independent scaling of the near- and mid-IR regions. Right: Optical depth r (at z = 0.1) versus 
gamma-ray energy in TeV for each EBL scenario tested. The optical depths for the baseline and 
two additional EBL models shown in the left panel are shown as well. 

EBL. 

The models tested here follow a more observational approach and were derived from what will 
be referred to throughout the text as the baseline model. This baseline model follows the shape 
outlined by lower limits derived from galaxy counts obtained with the Hubble Space Telescope 



(Gardner et al. 2000 Madau & Pozzetti 2000), the Spitzer Space Telescope (Fazio et al. 2004 



Papovich et al.|2004 ), and the Infrared Space Observatory ( Elbaz et al.|2002" ). An approach similar 



to that of Mazin Sz Raue (2007), using third order splines, was implemented to produce this baseline 



EBL model. The advantage to this technique is that it allows one to easily adjust the shape of the 
SED through manipulation of the control points defining the spline curve. 

Since the absorption of TeV gamma-rays by EBL photons is sensitive to the intensity at near- 
and mid-IR wavelengths, the regions between ~0.3— 7 ^m and ~7— 50 ^um were independently scaled 
to explore a variety of EBL shapesj^ The left panel of Figure [l] illustrates the range of scenarios 
investigated (shaded region). The baseline model, indicated by the thick solid curve, is shown 
along with two additional models (dashed and dotted curves) illustrating various levels of scaling 
in the two wavelength regimes. The right panel of Figure [T] shows the calculated optical depths r for 
gamma-rays given each EBL scenario. All optical depths were calculated assuming the cosmological 
parameters Hq = 70kms~^ Mpc~^, 0^ = 0.3, and = 0.7. 



^The near and mid-IR regions were scaled as a whole but future work will investigate scenarios using a finer 
resolution (e.g, including a near-IR excess due to Population 111 stars ( Dwek et al.|2005 l). 
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To cover the full range of intensities, 18 scaling factors were used in the near-IR regime and 
29 in the mid-IR. This resulted in a total of 519 EBL scenarios]^ The scaling factors were chosen 
to be linearly distributed in logio{i'Iu). 



3. Determining the EBL 



The analysis methods used here to constrain the EBL are described in Sections 3.1 and 3.2 In 
Section 3.3 we describe how these two methods were combined to derive even stronger constraints 
on the allowable EBL intensity parameter space. 



3.1. Method 1 - Spectral Shape Method 



With the advent of the Fermi Large Area Telescope (LAT) ( [Atwood et al. 2009), high energy 
observations of blazars are now possible in a regime where EBL attenuation is minimal. The 
operational energy range of the LAT spans ~0.1— 300GeV, overlapping observations with Imaging 
Atmospheric Cherenkov Telescopes (lACTs) ranging from ~0.05 — lOTeV (Hinton 2004 Albert 



et al. 2008 Holder et al. 2008). Since gamma-rays below ~10GeV travel unimpeded by the EBL 



(Stecker et al. 2006 Franceschini et al. 2008 Kneiske et al. 2004 Gilmore et al. 2009), the first 



method described here used spectra measured by Fermi as a proxy for the intrinsic spectra in the 
TeV regime. Measurements by lACTs were used to calculate a de-absorbed spectrum for each EBL 
model discussed in Section [2l 



The intrinsic TeV source spectra were determined using the relationship 



dEJ 



intr 



dEJ 



(1) 



obs 



where (dA^/d£')intr is the intrinsic spectrum, (dA^/d-E)obs is the observed spectrum, and t{E,z) is 
the optical depth at energy E and source redshift z. The intrinsic spectrum was assumed to have 
a power-law form given byj^ 

where No is the normalization at energy Eq, E is the energy, and T is the spectral index. The 
condition used to evaluate the consistency between the intrinsic TeV spectrum and the extrapolated 



■^The total number of scenarios would have been 522 if not for an imposed restriction that the mid-IR intensity 
be less than the near-IR intensity. This resulted in the exclusion of three models. 

power-law is generally a good approximation to blazar spectra over a limited energy range (e.g., 0.1 GeV < 
E < 100 GeV or 0.1 TeV < < 10 TeV) for measurements with limited statistics. 
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Fermi spectrum was given by 



|rT< 



GeV 



< N 



a, 



TeV 



GeV 



(3) 



where FxeV and cr^gv are the calculated lACT intrinsic spectral index and variance, respectively, 
TceV and CQgY are the Fermi spectral index and variance, respectively, and is the confidence 
level in units of standard deviations. An EBL scenario was considered viable if the de-absorbed 
spectrum satisfied the set criteria for consistency with the intrinsic spectrum predicted by the Fermi 
extrapolation. This criterion was that FxeV be within ±Na of FceV (i-e.. Equation |3]) , where N is 
1,2, or 3. 

As indicated by Equation [3j it was assumed that hard spectrum blazars have no intrinsic 
spectral break between the Fermi and lACT energy regimes up to several TeV. This assumption 
was motivated by both observational and theoretical evidence. Firstly, observations of nearby hard 
spectrum blazars, such as Mrk 421 and Mrk 501, indicate that their observed IC peaks can be 



located at energies >lTeV (Dwek & Krennrich 2005 Aleksic et al. 2010). Since EBL absorption 



softens the observed spectrum with respect the the intrinsic one, accounting for any absorption 



whatsoever will only move this peak to higher energies. In fact, Dwek Sz Krennrich (2005) show 



that most EBL realizations result in an intrinsic peak energy between ~1 and 5 TeV for the blazar 
H 1426+428 {z = 0.129). It is also well known that the spectra of Mrk 421 ([Krennrich et al. 



2001 Aharonian et al. 2002) and Mrk 501 (Samuelson et al. 1998: Aharonian et al. 1999) have 



exponential cutoffs at energies >4TeV, further supporting the assumption of an IC peak located 
at TeV energies. 

Population studies of blazars performed with Fermi show that the hardness of the Fermi 
spectrum is correlated with the synchrotron and IC peak frequencies (Abdo et al. 2010a|b ). In 
other words, a harder Fermi spectrum generally means a higher IC peak energy. The sources used 
in Method 1 have some of the hardest spectra of all Fermi-detected blazars (even harder than 
Mrk 421 and Mrk 501) and hence are likely to have some of the highest energy IC peaks. This 
provides an additional motivation for the assumption that, for these particular sources, the intrinsic 
spectrum up to a few TeV is consistent with an extrapolation of the spectrum measured by Fermi. 



Costamante et al. (2003), Katarzynski et al. (2006), Aharonian et al. (2007a), and Krennrich 



et al. (2008) have all shown that extremely hard intrinsic blazar spectra are unavoidable, for 



these already hard sources, using current EBL model estimates and observationally derived limits. 



Katarzynski et al. ( 2006 ) demonstrated that an IC peak at very high energies can be incorporated 



into the standard synchrotron self-Compton (SSC) emission model of blazars by introducing a 
low energy cutoff in the parent electron distribution responsible for the synchrotron and SSC 



emission. Furthermore, Tavecchio et al. (2009) showed that recent measurements from Swift provide 



supporting evidence for a very high energy IC peak in the hard spectrum blazar lES 0229+200. 



The spectrum of lES 0229+200 is the only one used here extending to energies that could 
potentially span the blazar IC peak for the case of hard spectrum blazars. As such, the analysis of 
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lES 0229+200 warranted a slightly different treatment |^ As shown in Figure 2(d) , the spectrum of 
lES 0229+200 ( [Aharonian et al.||2007d ) spans the energy range of ~0.5— 15TeV, nearly 1.5 orders 
of magnitude. This broad range, combined with the measurements above several TeV, necessitates 
the inclusion of a term to account for curvature in the intrinsic spectrum. This was done using a 
log-parabolic function of the form 



Ie) 



Nn 



intr 



E 

Eq 



-r-l3log^o(E/Eo) 



(4) 



where Nq is the normalization at energy Eq, E is the energy, T is the spectral index, and /3 is 
the log-parabolic fit parameter. Additionally, an F-test was performed, along the lines of Dwek &: 



Krennrich (2005), to test for the presence of an exponential rise in the intrinsic spectrum at the 



highest energies. This exponential fit had the form 



dEJ 



intr 



(-) 



(5) 



where {dN/dE)iag is the log-parabolic function from Equation |4| 7 is the exponential fit parameter 
(with units of energy), and E is the energy. If the resulting P- value from the F-test was > 95%, the 
exponential function was considered to yield a statistically significant improvement to the fit. The 
EBL model was then excluded based on the fact that it produced an unphysical intrinsic spectrum. 



A more detailed discussion of the F-test and its application is found in Dwek & Krennrich (2005) 



Using the analysis described above, the hard spectrum blazars lES 0229+200 (Aharonian et al. 



2007d), RGB J0710+591 (Acciari et al. 2010b), lES 1101-232 (Aharonian et al. 2006), and lES 



1218+304 (Acciari et al. 2010a[ ) have been studied. The combined Fermi and lACT spectra for 
these sources are shown in Figure [2j It is clear that the spectral breaks between the GeV and TeV 
regimes, as well as the energy ranges covered, vary between sources. This results in the derivation 



of a unique set of EBL constraints for each spectrum (Section 4.1 ). Hard spectrum blazars are ideal 



for these studies because they are more readily detected at TeV energies (Aharonian et al. 2007d 



Krennrich et al. 2008). They also stand the greatest chance of being detected at energies upwards 



of 10 TeV, broadening the EBL constraints into the mid-IR regime. Combining the analyses of 
multiple spectra further constrains the allowable parameter space. In the following, we refer to this 
approach as Method 1, or the spectral shape method. 



3.2. Method 2 - TeV Spectral Break Method 



The absorption of gamma-rays by EBL photons may produce breaks in the observed blazar 
spectrum for some SEDs. This break is produced by changes in the slope of the 7-7 optical 



*It should be noted that this "different" treatment was also applied to the other sources used in Method 1. 
However, it had no impact on the final results. 
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(a) lES 1218+304 (VERITAS ( jAcciari et al.|2010a| )). (b) IBS 1101-232 (H.E.S.S. ( [Aharonian eraL][2006l )). 
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(c) RGB J0710+591 (VERITAS ( |Acciari et al.|2010b| )). (d) lES 0229+200 (H.E.S.S. ( [Aharonian et al.|2007d| )) 



Fig. 2. — Combined Fermi and lACT E'^{dN/dE) spectra in units of GeVcm~^s~^. The best fit 
Fermi spectrum is indicated by the shaded region. The Fermi flux points were calculated by fixing 
the spectral index to the value obtained from the fit over the full energy range and then fitting the 
integral flux over each energy bin. The lACT spectra are given by the red line and flux points in 
each plot. Note: Fermi has a weak detection of lES 0229+200 (the current result with ~19 months 
of data is ~4(t). For this analysis we have fixed its spectral index in the Fermi regime to a value 
of 1.5. 
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depth. As illustrated in the right panel of Figure [T| the optical depth calculated for certain EBL 
scenarios is nearly flat in the energy range of ~1— several TeV resulting in an approximately energy 
independent absorption of gamma-rays. The observed spectral index in this energy range will be 
closer to the intrinsic value - producing a break in the spectrum at ~1 TeV. The magnitude of 



this break increases with the source redshift and depends on the near- to mid-IR ratio (Imran & 



Krennrich 2008). With currently available data, no individual blazar spectrum is likely to show 
a statistically significant break in its spectrum. However, a large sample of blazars may reveal a 
redshift dependent trend. The presence/absence of such a trend in observations can be used to 
constrain the EBL. Table [T] lists the sample of blazars used to perform this study. 

The shape of the EBL uniquely determines the spectral break versus redshift distribution and 
was calculated here for a multitude of scenarios using a test blazar spectrum. This test spectrum 
was chosen to represent the "average" spectrum of the blazars listed in Table [T] This was motivated, 
for example, by the fact that the spectrum of lES 0229-1-200 ( Aharonian et alT]|2007d ) ranges from 



0.5 -15 TeV whereas that of lES 1218-^304 (Acciari et al. 2010a) spans ~0.2-2TeV. Using a 



test blazar spectrum ranging from 0.2— 15 TeV would be unrealistic as there is no individual source 
in the sample covering such a wide energy range. To generate a spectrum representative of the 
overall sample of blazars, the average lowest energy bin, highest energy bin, and number of bins 
in the sample were calculated. This resulted in a test spectrum characterized by 8 energy bins 
ranging from 200 GeV to 5 TeV. The flux normalization and spectral index (intrinsic values) of the 
test blazar spectrum are irrelevant since only the magnitude of the spectral break is of concern. 
Nevertheless, values need to be chosen to perform the analysis. The Crab Nebula flux at 1 TeV was 
used for the normalization along with an intrinsic spectral index of 1.5. The underlying assumption 
in this method is that the intrinsic blazar spectrum is well described by a single power-law over the 
energy range considered. 

Beginning with an exact power-law function, and calculating the spectrum due to EBL absorp- 
tion, one is left with an observed spectrum that shows detailed features of the 7-7 optical depth 
and is decidedly not of a power-law (or broken power-law) form. This is, of course, not what is 
observed by lACTs. To produce an absorbed spectrum more consistent with observations, each flux 
point was assigned an error of 25% of the flux in that bin. This value was chosen to approximately 
coincide with the mean error of all flux points from the lACT measurements in the blazar sample. 
Gaussian fluctuations were added to the spectral points using a Normal distribution with a stan- 
dard deviation equal to the 25% error bars. Optical depths for the test spectrum were calculated 
for all EBL scenarios over the redshift range 0.05-0.45 in steps of 0.05. No EBL evolution with 
redshift was implemented. 

The absorbed spectra were calculated using the inverse relation of Equation [T] Each spectrum 
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was fit with a broken power-law of the form 



dE 



^0 { ] , E < inbreak 

-^break / 



(6) 

No I , E > -Ebreak 

k V -C'brcak / 



where A^o is the normalization at the break energy inbreak , Ti and r2 are the spectral indices below 
and above inbreak) respectively, and E is the energy. The spectral break was defined as 

Ar = Ti - Ts . (7) 

With this definition, a spectral hardening above the break energy would result in a positive AF. 
From here the distribution of AF as a function of redshift was determined for each EBL scenario. 
This process was repeated 100 times, mimicking multiple independent observations, each time 
yielding a different result due to the Gaussian fluctuations in the spectra. This facilitated the 
determination of the break versus redshift distribution to an arbitrary precision and accuracy. 

Each distribution was fit with a line, i.e., 

Ar{z) = mz + b , (8) 

where AF(2;) is the spectral break at redshift z and m and b are free parameters. Using the best 
linear fit obtained from observational data, the la, 2a, and 3a two dimensional contours in slope 
and intercept space (i.e., m and b, respectively, in Equation [s]) were determined. The best fit 
results from the expected spectral break versus redshift distribution for each EBL scenario were 
then compared with the observationally derived contours. If the fit parameters for a particular 
scenario fell within the two dimensional la contour, this scenario was included as part of the la 
contour in EBL parameter space. The 2a and 3a contours in EBL parameter space were determined 



in the same fashion (Section 4.2). In the following, we refer to this approach as Method 2, or the 
TeV spectral break method. 



3.3. Complementarity of Methods 1 & 2 in Constraining the EBL 

The constraints placed on the EBL are presented in the form of a contour plot. This contour is 
drawn in the parameter space defined by the near-IR intensity at 1.6 //m and the mid-IR intensity 
at 15 /xm. These two wavelengths are representative of the positions of the near-IR peak and 
mid-IR trough in the EBL SED. Given that the EBL models tested are relatively smooth, these 
two parameters serve as an adequate characterization of a given EBL scenario. An illustration of 
such a contour plot is shown in Figure [3j Each point represents one of the EBL scenarios tested. 
Shown in this fashion, the linear sampling of intensity in log-space is clearly seen. A constant 
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(a) (b) 

Fig. 3. — EBL intensity at 15 fim plotted versus the intensity at 1.6 ;um. Each grey dot represents 
an EBL scenario tested, (a) As the intensities at 1.6 and 15 /xm increase (moving up and to the 
right on the plot), the inferred intrinsic blazar spectrum becomes harder. The intrinsic spectrum 
softens as the two intensities decrease (down and to the left). This is indicated by the arrows in the 
figure. By placing restrictions on the maximum hardness and softness of the intrinsic spectrum, one 
obtains a contour in intensity space oriented approximately perpendicular to these two directions 
(shaded region), (b) Increasing the near-IR to mid-IR ratio increases the redshift dependence of 
the EBL-induced spectral break at TeV energies (i.e., AT(z)/z becomes larger). Decreasing this 
ratio results in an increasingly negative dependence of spectral break versus redshift (i.e., AT(z)/z 
becomes smaller to the point of being negative). This is indicated by the arrows in the figure. By 
comparing the redshift dependence of the TeV spectral break produced by various EBL models 
with observations, one obtains a contour in intensity space oriented as shown by the shaded region. 



ratio of intensities (i.e., i^/,y(1.6^m)/z^/;y(15;um)) would be given by a straight line with positive 
slope. It should be noted that parametrizing the results in this fashion does not mean the derived 
constraints only apply to the EBL intensity at 1.6 /xm and 15 /im. This presentation is simply a 
way of characterizing the average slope of the EBL SED between these two wavelengths. 



The shaded region in Figure 3(a) shows how a contour, derived from the analysis outlined in 
Section 3.1, might appear. As the near- and/or mid-IR intensities increase (decrease), the inferred 
intrinsic spectrum hardens (softens). A contour oriented as in Figure 3(a) is obtained when placing 
limitations on the maximum hardness/softness of the intrinsic spectrum. 



An approximately orthogonal contour is obtained from the analysis outlined in Section 3.2 



This is because, when increasing the near-IR to mid-IR intensity ratio, the dependence of the 
EBL induced spectral break with redshift becomes increasingly positive. That is to say, AT{z)/z 
becomes larger. Conversely, decreasing the near-IR to mid-IR ratio makes AT(z)/z smaller to 
the point where it can become negative. These directions are shown on the contour illustration 
in Figure |3(b)[ By comparing the redshift dependence of the TeV spectral break produced by 
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different EBL models with that from observations, one can obtain a contour in EBL intensity space 
as represented in Figure [3(b) [ by the shaded region. 

Comparing the two contours from Figure[3| it becomes clear that the EBL constraining methods 
outlined in Sections 3.1 and 3.2 are complementary to one another. Each method is sensitive to 
different characteristics of the EBL SED thereby improving the constraints at near- and mid-IR 
wavelengths. 



4. Results 



The techniques described in Sections |3.1| and 3.2 have been applied to a sample of blazars 
(Table [T]) . The four blazars selected for the spectral shape analysis were chosen for their redshifts 
{z > 0.1), hard spectra, and steady emission in the Fermi regime. 

The blazar sample used in the TeV spectral break analysis was selected based on characteristics 
of the source spectra in the TeV regime. Most importantly, to test for a spectral break at ~1 TeV, 
each spectrum must have at least two data points both above and below the break energy. Spectra 
from large flaring events were excluded as their atypically high statistics skew the overall results. 
An example of this was the ~7Crab flare event of PKS 2155-304 (Aharonian et al. 2007c). The 



two nearby blazars Mrk 421 and Mrk 501 were also excluded from this analysis due to the fact that 
the presence of a cutoff at ~4TeV in their intrinsic spectra ( Krennrich et al.|2001 Samuelson et al. 
1998), combined with their low redshift (z ~ 0.03), greatly inhibits the search for a spectral break 
resulting from EBL absorption. 



The analysis results for Methods 1 and 2 are presented in Sections 4.1 and 4.2, respectively, 
and are then combined in Section [4.31 



4.1. Spectral Shape Analysis 

The source spectra used for the spectral shape analysis are shown in Figure [2] Each Fermi 
spectrum consists of approximately 19 months of data and was analyzed using the Science Tools 
v9rl5p2 package distributed by the Fermi collaboration|^ Only events with zenith angles less 
than 105° were included to avoid contamination from Earth's gamma-ray albedo. All exposures 
and fluxes were calculated using the post-launch instrument response function P6_V3_DIFFUSE. The 
galactic and extragalactic diffuse gamma-ray backgrounds were modeled using data available on 
the Fermi Science Support Center websitelj The variability index of each source, as quoted in the 



^http:/ /fermi. gsfc.nasa.gov/ssc/data/analysis/software/ 

®http:/ /fermi. gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels. html 
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Fermi 1 Year Catalog^ is less than 23.21 indicating a >99% probability the source emission is 
steady. Spectral fitting was performed from 200 MeV to the highest photon energy associated with 
the given source. TeV spectra were obtained from the references listed in Table [l| 

The Fermi analysis of lES 0229+200 was handled differently than for the other three sources 
due to the fact that Fermi has only a weak detection of this source. Nevertheless, it was included 
since it has the hardest TeV spectrum of all the blazars studied here and provides a unique set of 
constraints. To obtain a spectrum from the Fermi data, an assumed spectral index of 1.5 it 0.2 
was used. This parameter was fixed in the source model so that the the best fit flux normalization 
could be determined. The overall test statistic (TS) for this source, between 1 and 28GeV, was 
18.0 (~4(t). The TS values for the first (1—5 GeV) and second (5—28 GeV) energy bins were 2.6 and 
15.4, respectively. 

The la, 2a, and 3a contours for each blazar, constraining the EBL intensity at 1.6 ^m and 
15 /xm are shown in Figure |4] The reader is cautioned not to view these contours as strict exclusion 
regions for other authors' models since the analysis is dependent on the shape of the EBL SED. 
However, given that most models have the same general shape, it is a safe assumption that scenarios 
falling outside the 3a contours derived here are not consistent with this analysis. 

Figure |4] shows that gamma-ray observations provide strong constraints on the 1.6 /im EBL 
intensity, but leave the 15 fini intensity largely undetermined. This is due to the lack of TeV data 
extending up to multi-TeV energies which provide the strongest constraints on the mid-IR region 
of the EBL. It is worthwhile to point out some particular properties of the different contours in 
Figure [4] as they provide insight into the constraints derived with each blazar. 

One of the most noticeable features is the narrowness of the lES 1218+304 contours (Figure 
Qa)). This is a result of the high precision measurement of the spectral index by Fermi and 
VERITAS in both the GeV and TeV regimes (see Table [T]) . The measurement errors are a factor of 
2 or more smaller than those of the other blazars in the sample. Another distinguishing feature is 
the vertical orientation of the contours indicating that the TeV spectrum is minimally affected by 
EBL absorption in the mid-IR regime. This results from the fact that the spectrum only extends 
up to ~2 TeV. 

The spectrum of lES 1101-232 extends up to ~3TeV but is also insensitive to the mid-IR 
portion of the EBL (Figure |4]^b)). This is most likely due to the somewhat large error on the 
spectral index measured by Fermi (1.87 ± 0.28) as well as the large power-law fit probability 
(~80%) to the observed TeV spectrum (see [Aharonian et al.| (|2006[))f1 



^http:/ /fermi. gsfc.nasa.gov/ssc/data/access/lat/lyr .catalog/ 

® A fit probability of 80% suggests the possibility of either overestimated errors or somewhat correlated data points. 
In either case, fluctuations in the highest energy spectral points (those sensitive to the mid-IR portion of the EBL) 
will have a limited impact on the overall power-law fit to the de-absorbed spectrum. 
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Fig. 4. — Constraints, provided by Method 1, on the EBL intensity, in units of nW m^^ sr^-*^, at 
1.6 fim and 15 ^m. The contours shown are for the 1 (red), 2 (blue), and 3 (green) sigma confidence 
intervals. Each grey dot represents one of the EBL scenarios tested. Note: Fermi has not obtained 
a significant detection of lES 0229+200. For this analysis we have assumed the spectral index in 
the Fermi regime to be 1.5 ± 0.2. 
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In contrast to lES 1218+304 and lES 1101-232, the contours for RGB J0710+591 (Figure 
|4]^c)) are somewhat incUned, indicating sensitivity to the mid-IR EBL intensity. The spectrum of 
RGB J0710+591 extends up to ~4TeV demonstrating that a modest increase in the maximum 
energy of the spectrum can have a significant impact on EBL studies. 

The spectrum of lES 0229+200 is unique in that it starts at ~500 GeV and extends up to 
~15 TeV. As a result, this source is the most sensitive, of ah those analyzed here, to absorption by 
the mid-IR portion of the EBL. This is evident in Figure Qd) . While none of the aforementioned 
blazars strongly constrain the 15 /um EBL intensity, the spectral shape analysis of lES 0229+200 
constrains the intensity to be <5 nWm~^ sr~-^. 



4.2. TeV Spectral Break Analysis 



The energy dependence of the 7-7 optical depth for some EBL scenarios may produce a break 
around ~1 TeV in blazar spectra (Section 3.2 ). Here, a spectral break energy of 1.3 TeV is assumed. 
The choice of this value is easily arrived at using the the fact that the 7—7 cross section peaks at 
Xe{iJ,iii) ~ 1.24£'^(TeV), where is the EBL photon wavelength and is the gamma-ray energy 
( |Dwek Sz Krennrich 2005| ). Since the near-IR peak of the EBL occurs at ~ 1.6 fim, a flattening 
of the 7-7 optical depth would be expected at 1.3 TeV. The right panel of Figure [T] shows 

that the optical depth does, in fact, flatten out between 1 — 2 TeV, further motivating this choice 
for the spectral break energy. 

To determine the spectral break versus redshift distribution {AT{z)), a sample of 12 blazars 
(Table [T]) with spectra extending both above and below the break energy were fit with a broken 
power-law (Equation [g]) . Figure [s] shows the distribution of Ar(z) versus the source redshift for 
all blazars in the sample. Two linear fits were performed on the data - one leaving the intercept 
Ar(0) as a free parameter and one with the intercept fixed to 0. The results of these two fits were 



AT{z) 



' (8.68 ± 5.37)z - (0.24 ± 0.71) , y^jv = 14.10/10 
(6.95±1.65)z , xV'^ = 14-22/11 



(9) 



In the absence of any spectral break one would expect Ar(2;) 
at the ~3.2(T level. 



0. The data exclude this hypothesis 



Figure l5l also shows the Ar(2;) distributions produced by the models of Aharonian et al. 



r scaled by a factor of 0.55, and Franceschini et al. (2008). A summary of how well these 



models describe the data, along with the other authors' models considered here, is given in Table 
[2j Column 3 gives the best linear fit to the Ar(2:) distribution produced by the given model. This 
is derived in the same way outlined in Section |3.2[ However, in this case we use the EBL model 



^The EBL SEDs from 



Aharonian et al. 



( 2006 1 are scaled versions of the model from 



Primack et al. 



(20011 
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from the reference listed rather than one of our own reahzations discussed in Section [2j With the 



exception of the EBL scenario from (Aharonian et al. 2006), all models are inconsistent with TeV 
spectral break data at the >3a level. 

Figure |6] shows the la, 2a, and 3a contours resulting from comparison between the Ar(z) dis- 
tributions generated from the EBL models tested and observations (Figure [s]) . In addition, points 



taken from Aharonian et al. (2006), Kneiske et al. (2004), Stecker et al. (2006), Franceschini et al 



generated here 
tion 



10 



(2008), and Dommguez et al. (2011) show where these other EBL models fall within the contours 



The orientation of the contours are as anticipated from the discussion in Sec- 



3.3 The limits on the near- to mid-IR ratio (z^/jy(1.6 ^m)/z^/jy(15 /xm)) become increasingly 



restrictive as the intensities at 1.6 /xm and 15 fim increase. Using this technique, future direct mea- 
surements of the EBL at near-IR wavelengths (for instance by the James Webb Space Telescop^^ 
will help to further constrain the EBL in the mid-IR regime. 



4.3. Combining the Analyses 

Improved constraints on the EBL can be obtained by combining the results of the analyses 
from Sections 4.1 and 4.2 Figure [7] shows the overlay of 2a (a) and 3a (b) contours from the 
analysis of RGB J0710+591, lES 1218+304, and Ar(z). The contours from lES 1101-232 and 
lES 0229+200 are left off because they do not further restrict the constraints. Figure [8] shows the 
final allowed 2a and 3a regions of parameter space. The overlapping region indicates a very low 
EBL intensity at 15 fiui. It should be reiterated that the assumption here was that the intrinsic 
TeV spectrum was approximately equal to the Fermi spectrum. If, in fact, there is an intrinsic 
break between the GeV and TeV spectra of blazars, the intrinsic TeV spectrum will be softer than 
is assumed here. This will move the Fermi-IACT derived EBL contours toward the left of Figure 
[sjand hence allow for lower EBL intensities at 1.6 /^m. However, given the Ar(z) derived contour 
(Figure [6]), the EBL intensity at 15 /im will then be constrained to even lower values. 



Figure [9] shows the 2a confidence region for the EBL SED. The models of Aharonian et al. 



(2006), Stecker et al. (2006), Franceschini et al. (2008), and Kneiske et al. (2004) are also shown. 



This region is generated using the rectangular area encompassing the 2a contour in Figure ^ 



It should be noted that some of these EBL models have bumps at ~15 fj,m. The data points should therefore only 
be used as a point of reference and not for a strict comparison. The compatibility of these models with our results 
was determined using the full EBL SED. 



'http:/ /www .jwst.nasa.gov/ 
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5. Caveats and Considerations Regarding Methods 1 & 2 

There are a number of caveats pertinent to the analyses presented here. While each of these 
cannot be discussed exhaustively, we address the most relevant topics below. First, however, we 
make a few general remarks. 

The EBL attenuation factors calculated throughout the paper were done so using the central 
energy (in logarithmic space) of each spectral energy bin. Another way of doing this would be 
to calculate the average attenuation across the entire energy bin. We have verified numerically 
that these two approaches result in absorption factors within ~1% of one another given the typical 
energy bin widths encountered here. 

Not every possible shape of the SED was explored with the EBL parametrization used here. 
The aim of this work is to constrain the overall shape of the EBL and not to test for various bumps 
and wiggles. The reader is therefore cautioned not to view any individual data point, derived from 
galaxy counts, in Figure [9] as being inconsistent with our results if it lies outside the designated 
contour region. This analysis does not rule out the possibility that there are, for instance, bumps 
in the EBL spectrum (e.g., at ~15 ^um). What is demonstrated by this analysis is that the near- to 
mid-IR ratio of the EBL appears to be larger than is predicted by many models thereby requiring 
a relatively high intensity in the near-IR (see Section u|. 



5.1. Hard Spectrum Blazars and TeV IC Peak Energies 

The assumption made in Method 1 was that hard spectrum blazars have IC emission com- 



ponents peaking at TeV energies. We believe this assumption is well motivated (see Section 3.1) 
but, if it were to prove incorrect, and one instead allowed for some level of intrinsic curvature, the 
lower limits obtained here in the near-IR would be pushed to lower intensities. One could also 
interpret the upper boundary of the 2a contour in Figure [9] as an upper limit. However, there is a 
very important caveat to this second interpretation. If the near-IR intensity decreases, the mid-IR 
intensity must also decrease so that the overall slope of the EBL SED between near- and mid-IR 
wavelengths remains within the allowed range determined with Method 2. 

It should also be noted that the la contours shown in Figure |4] do not all overlap in a common 
region of EBL parameter space. This could be interpreted as being due to intrinsic differences 
between some of the blazars (e.g., different intrinsic spectral breaks). However, given that this is 
a la effect, and a common overlapping region for the 2a contours can be identified, we do not 
consider this issue to be of significance. 
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5.2. Spectral Break Versus Redshift Calculation 



There are multiple ways of to go about calculating the spectral break versus redshift distribu- 



tion for a particular EBL scenario. The method implemented here (Section 3.2) used a test blazar 



spectrum representative of the average spectrum of the 12 blazar data sample. This determination 
governed both the energy range covered and the errors on each flux point. An absorbed spectrum 
was then calculated for a uniformly spaced set of redshifts. 

An alternative to this approach is to use 12 different test blazar spectra, with each one being 
representative of a particular blazar in the dataset. In this case, one would calculate the absorbed 
spectrum, at each redshift of the observed blazars, using a test spectrum with energy bins and error 
bars equivalent to the associated observed blazar spectrum. This technique is most appropriate 
when testing EBL models that contain features in the SED (e.g., the bump at ~15//m in the model 



of Dommguez et al. (2011)). In cases such as this, the spectral break calculated, due to EBL 



absorption, can be sensitive to the energy range over which one performs the fit. Using the same 
energy range and binning of the spectral observations allows one to determine the expected spectral 
break given the specific dataset and an assumed EBL scenario. This technique is also preferable 



when the mid-IR EBL intensity is quite high (e.g., the Stecker et al. (2006) fast evolution model) 



We have used the technique described above to compare the various EBL realizations of other 
authors with the spectral break versus redshift data. The exclusion probabilities obtained are 
shown in Table [2] along side those calculated using the simple best linear fit to the spectral break 
versus redshift distribution for each scenario. With the exception of the fast evolution model of 



Stecker et al. (2006), the values obtained using the two techniques are comparable. The reason for 



the 1.5(T difference in the results of the Stecker et al. (2006) fast evolution model is likely due to the 



fact that it has the highest EBL intensity, and therefore the largest amount of absorption, which 
makes it the most sensitive to the energy range over which the absorbed spectra are fit. 

As an additional cross check of Method 2, we performed the exact same analysis described in 



Section 3.2 using a test blazar spectrum with 6 spectral points ranging in energy from 400 GeV 



to 3 TeV. This more limited energy range around 1 TeV elucidates how the low and high energy 
portions of the spectrum affect the spectral break calculation. The contours obtained using this 



test spectrum are shown in Figure 10 The main effect of the reduced energy range is to shift the 
contours upward to slightly higher mid-IR intensities. This is due to the fact that, by excluding 
energies above 3 TeV, the spectral fit is less impacted by the sharp rise in absorption, above a few 
TeV, that can occur for high mid-IR EBL intensities. This results in the calculation of a spectral 
break versus redshift distribution with a more positive, or less negative, slope which allows for 
higher mid-IR intensities. The corresponding 2a contour region for the EBL SED is shown in 
Figure [9] as the partially transparent blue shaded region. 

It is also important to understand the effects of each individual source on the observed spectral 
break versus redshift distribution shown in Figure [5] Namely, how do the results change if sources 
lying >2(T outside the fit (e.g., PKS 0548-322 and H 1426+428) are removed? By excluding H 
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1426+428 from the dataset, the fit result of Figure [s] changes by only ~5%. Removing PKS 0548- 
322 from the dataset results in a best fit of r{z) = (12.68 ± 5.79)z - (0.85 ± 0.78). In this case, 
both the slope and intercept of the fit change by less than la while the exclusion significances of 
the various EBL models mentioned in the paper (e.g., Primack et al. (2005), Stecker et al. (2006), 
etc.) change only marginally. 



5.3. Evolution of the EBL with Redshift 



The evolution of the EBL as a function of redshift should be taken into account when calculating 
the gamma-ray absorption for sources with large redshifts. Since the quality of the data used, and 
the nature of the analysis performed, will impact exactly at what redshift EBL evolution becomes 
a necessary inclusion, we investigate its effects here. 

The EBL is comprised of photons generated by sources that have been evolving since the early 
Universe. This evolution is, therefore, imprinted on the EBL but has been neglected in our analysis 
thus far. This means that, when calculating the absorption for a source at a redshift of, say, z = 0.2, 
we have been assuming that all of the target EBL photons are in place at this redshift and that 
there is no absorption and/or reemission of these photons (i.e., no evolution) between z = 0.2 and 
z = 0. This assumption is reasonable for small redshifts, but one must determine at what redshift 
it is no longer warranted. 

The evolution of the EBL has been examined in detail by many of the authors whose models 
are used here (e.g., Kneiske et al. (2004) and Franceschini et al. ( 2008| )). There is also a simple 
first order approach to incorporating EBL evolution into absorption calculations. This method 



simply applies a correction factor to the cosmological photon numbers density scaling (e.g., Madau 



Sz Phinney (1996) and Raue Sz Mazin (2008)). This factor is chosen such that the resultant EBL 



evolution is in good agreement with more detailed models (e.g., Kneiske et al. (2002)). [Raue &: 



Mazin (2008) have shown that this approach can be used to match the EBL evolution well for 



redshifts of z < 0.7. 

Using this approach, we have investigated the effects of EBL evolution on the results presented 
here. The inclusion of EBL evolution lowers the calculated optical depths by ~10% at a redshift 
of z = 0.2. For the blazars analyzed in Method 1, this translates to a softening of the intrinsic 
spectra by ~0.5(T. This corresponds to a shift of ~10%, to higher EBL intensities, in the contours 
of Figure |4j 

Method 2 is insensitive to the intrinsic spectral softening effects due to the incorporation of 
EBL evolution since this method relies on the measurement of the spectral break and not the 
spectral index itself. In other words, it depends on a relative measurement rather than an absolute 
one. It should be noted, however, that this is contingent upon the (reasonable) assumption that 
only the intensity of the EBL evolves over such small redshifts and not the shape of the SED. 
In other words, the relative contributions to the EBL from different emission components remain 
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approximately the same over small redshifts. 



5.4. Alternatives Scenarios Describing Blazar Spectral Observations 



Alternative scenarios that could explain the observed spectra of blazars, while allowing for a 
wider range of EBL intensities, have been investigated by others. Some have proposed that the 
observed high energy gamma-ray emission from blazars could be dominated by secondary photon 



emission produced along the line of site by cosmic-ray interactions with EBL photons (Essey & 



Kusenko 


2010 


Essey et al. 


2010a|b 



used to describe the observed spectra of many VHE blazars. Other authors have suggested that 
photons could oscillate into light axion-like particles (ALPs) via interactions with intergalactic 



magnetic fields (De Angelis et al. 2007, 2009). These ALPs, which travel unimpeded through 
the Universe, can then convert back into photons before reaching the Earth. Since ALPs are not 
affected by EBL absorption, this process would allow for an EBL intensity that is, in actuality, 
higher than would appear based on standard absorption calculations. These are just two of the 
many interesting avenues of research being pursued in relation to the EBL. 



6. Improving Constraints on the EBL 

The constraints on the EBL presented here can be improved through further observations 
with Fermi and lACTs. Fermi is continuously monitoring the entire sky and hence data on the 
sources used in Method 1 are being acquired all the time. Improved measurements in the TeV 
regime using current instruments will require either deep or flaring state observations of blazars 
and will predominantly impact the constraints obtained with Method 2. As can be seen from 
Figure [5j the majority of the blazar spectral breaks measured using current data are not, on an 
individual basis, statistically significant. The driving factor behind the large measurement errors is 
an insufficient knowledge of blazar emission at energies >1 TeV for most of the distant sources. As 
deeper exposures are obtained, and source photon statistics continue to improve, measurements of 
many blazars will extend to ever higher energies until their spectra either intrinsically, or via EBL 
absorption, cut off. 

As spectra span a larger energy range, and as the precision of each flux point improves, it is 
less likely that their intrinsic emission will be well described by a single power-law. When looking 
for evidence of a break at ~lTeV, it would be necessary to modify the approach used here by 
either truncating the fit regime so as to exclude energies above or below where the IC peak is 
believed to be located, or use a more complicated fit function that includes (at least) a spectral 
break parameter along with a curvature parameter. It will also likely be necessary to use the more 



sophisticated spectral break versus redshift determination described in Section 5.2 



Another consideration, regarding the detection of blazars in flaring states, is that of spectral 
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variability. Some blazars exhibit spectral hardening with increasing flux (Krennrich et al. 2002 



Albert et al. 2007b|c ) while others show no evidence for spectral variability when flaring (Acciari 



et al. 2010a). Any change in the overall slope of the spectrum would have a limited impact on 



Method 2, given that it relies on the relative change in spectral index (spectral break), but the 
limits obtained using the Fermi and lACT data as in Method 1 could be affectedj^ If the TeV 
spectrum were to harden, while the Fermi spectrum remained unchanged, tighter constraints on the 
EBL could be obtained due to the smaller observed spectral break. However, if spectral variability 
takes place entirely above or below the IC peak, for hard spectrum blazars, where we have assumed 
that the IC peak is at multi-TeV energies, both Fermi and lACTs should see the same amount of 
spectral hardening thereby having no impact on the derived constraints. Simultaneous Fermi and 
lACT data are needed during both flaring and quiescent states to truly evaluate the effect (if it 
exists) of spectral variability on hard spectrum blazars. 



7. Discussion & Conclusions 

We have presented constraints on the EBL intensity using two distinct methods for the analysis 
of blazar spectra. The spectral shape method (Method 1) used Fermi observations of the three hard 
spectrum blazars RGB J0710+591, lES 1101-232, and lES 1218+304 as a proxy for the intrinsic 
source emission in the sub- TeV to TeV energy regime. In addition, lES 0229+200, which has a 
weak Fermi detection (~4cj), was used in this analysis, with an assumed Fermi spectral index of 
1.5, as its elucidates the impact of sources detected above ~10TeV on EBL constraints. Observed 
lACT spectra were corrected for attenuation due to different EBL realizations. Limits on the EBL 
were then placed by rejecting realizations that did not produce deabsorbed (i.e., intrinsic) lACT 
spectra consistent with the Fermi extrapolation. 

The second method, the TeV spectral break method, involved two steps. The first consisted of 
characterizing the spectral break, occurring at ~1 TeV, using the lACT observations of 12 blazars 
located at different redshifts. The magnitude of the spectral break was then plotted as a function 
of redshift (Figure [5]). We then assumed that the intrinsic blazar spectrum can be characterized as 
a power-law which, unlike Method 1, was not determined by the Fermi spectrum, and is therefore 
more general in nature. The second step consisted of comparing the data in Figure [5] with the 
redshift distribution of spectral breaks derived from the attenuation of an intrinsic test spectrum 
using different EBL realizations. 

We showed that Method 1 primarily constrains the overall EBL intensity, whereas Method 2 
primarily constrains the relative EBL intensities in the near- and mid-IR (~1.6/xm and ~15//m, 
respectively). The two approaches are therefore complementary and provide a closed confidence 
range for the EBL intensity. This is the first time these methods have been used to derive such a 



The four sources used in Method 1 show no signs of spectral variabihty. 
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contour for the EBL SED. 



The resulting constraints imposed by both methods on the the 1.6 /xm and 15 fim intensity of 
the EBL are shown in Figure l8| Using ultra deep 15 /um mapping observations of the gravitational 



lensing cluster Abell 2218 with the AKARI infrared satellite, Hopwood et al. (2010) derived an 



intensity of 1.9 it 0.5nWm~ sr~ for the integrated galaxy light (IGL) down to an intensity of 
~0.01 mJy. A model fit to the differential contribution to the IGL, combining data from several 
fields, gives a total integrated EBL intensity of 2.0 it 0.4nWm~^sr" -1 dHopwood et"aI1|2010[ ). The 
IGL seems to have totally resolved the EBL at 15 /um. Using our constraints, and the claimed 
15 /xm EBL detection, we derive a value of 17 it 3 nW m~^ sr~^ for the EBL intensity at 1.6 //m, the 



error representing a la uncertainty. Comparison with the IGL (Madau & Pozzetti 2000) suggests 
that ~50%, and perhaps as much as ~90% (Keenan et al. 2010), of the EBL has been resolved at 
this wavelength. 

Figure |9] presents the 2a confidence range of the EBL intensity derived in this paper as a 
function of wavelength (shaded blue). Our derived constraints are consistent with claimed EBL 
detections at 3.5 //m dPwek k Arendt||1998D and 3.6 //m ( [Levenson Wright||2008l ) . Figure |9] also 
shows that our EBL constraints are consistent with some models, such as the scaled version of 



Aharonian et al. (2006) and the most recent model of Hopwood et al. (2010). However, the models 



of Stecker et al. (2006), Franceschini et al. (2008), Kneiske et al. (2004), and Dommguez et al 



(2011) are disfavored by our analysis at >3<t confidence level. 



By inspection of Figure[9j one may be tempted to conclude that the models of [Franceschini et al. 



(2008) and Dommguez et al. (2011) are inconsistent with the analysis presented here predominantly 



because their intensities in the near-IR are too low. One might then suspect that this is purely due 
to constraints placed on the EBL by Method 1, which would be lower if the assumption of a single 
power-law between GeV and TeV energies was not made. This, however, would be an incorrect 
conclusion. The exclusion probabilities calculated in Table [2] are entirely independent of Method 
1, using only the spectral break versus redshift distribution of Method 2. As such, what is truly 
driving the agreement or disagreement of a particular model with the data in this case is the near- 
to mid-IR ratio. 

The minimum average slope between 2 and 10 fim, falling within our 2a contour in Figure |9j 
is a ~ 1.36 (ocA~°). This value is in agreement with the limit placed by Aharonian et al. (2007d) of 



a > 1.1 ± 0.25 using lES 0229+200 data from H.E.S.S. [Aharonian etaL] ( |2007dp also calculate an 
upper limit on the EBL intensity at 10 /xm of 3.1 nW m~^ sr~^. This constraint is also compatible 
with our results, which places an upper bound at 10 ^m of ~2nWm~^ sr~^. 

In summary, our analysis has shown the advantages of using these two complementary ap- 
proaches in constraining the EBL. It has also provided evidence that the near- to mid-IR ratio of 
the EBL may be larger than previously thought. Future measurements with Fermi, as well as both 
current and next generation lACTs, such as the Cherenkov Telescope Array (CTA), will continue 
to improve these constraints. 
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A. Investigation of Systematics 

To obtain a more thorough understanding of how the choice of certain parameter values (flux 
point error bars, number of iterations, redshift binning, etc.) affect this analysis, we have investi- 
gated the potential systematic effects introduced by each. 



A.l. Test Blazar Spectrum Flux Point Errors 



The application of 25% error bars to the test blazar spectrum (Section 3.2), while motivated 
by the mean percent error from VERITAS spectra, is a somewhat arbitrary choice. However, the 
systematic affect on the final result is limited, as illustrated in Figure [TT} This figure shows the fit 
results for the TeV spectral break versus redshift distribution (Equation [s]) , after 100 iterations, 
assuming a particular EBL model. The left (right) panel shows the best fit slope (intercept) using 
six different size error bars, ranging from 15%"65%, for the test blazar spectrum. While the error 



in the best fit value increases with the error in the test blazar spectrum, the standard deviation 
in the final result is <3% of the error in the best fit determination from observations (Figure [5]). 
From this one can be conclude that the systematic effect introduced by the choice of error bars in 
the test spectrum is minimal. 



A. 2. Redshift Binning and the Number of Iterations for AT{z) Determination 

The choice of binning in redshift and the number of iterations used for the determination 
of AT{z) are two potential sources of systematics. A redshift binning that is too course may 
result in a poor characterization of the spectral break versus redshift distribution. Many of the 
blazars used in the TeV spectral break analysis are separated in redshift from another source in 
the sample by Az < 0.05 (see Table [T]). It may seem then that a redshift binning of 0.05 is 
insufficient for comparing the AT{z) distributions for given EBL models with that from observation. 
Additionally, the choice of 100 iterations (i.e., simulated observations) for the determination of the 
Ar(z) distribution is not necessarily sufficient for the fit to converge to its true value. 



Figure 12 addresses both these points. The four panels show the spectral break versus redshift 



This preprint was prepared witii tiie AAS lAT^njX macros v5.2. 
^^This could be remedied by using a larger number of iterations (i.e., simulated observations) 
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distribution after 1, 25, 50, and 100 iterations using one particular EBL model. It can clearly be 
seen that the distribution has converged after 100 iterations. This is reflected not only in the data 
point error bars, but in the errors on the best fit parameters as well. It is also clear that, while the 
distribution is not perfectly linear, there is no evidence of features from small variations in redshift. 
In other words, a finer binning in redshift would have virtually no impact on the fit results. 



A. 3. Potential Pileup Effect at Very High Energies 

The possibility exists that a redshift dependent spectral break at ~1 TeV could be the result 
of an effect intrinsic to the observing instrument. Perhaps the highest energy spectral bins contain 
a pileup of events that systematically hardens the spectrum. We have investigated this issue by 
examining two distributions. The first distribution of relevance is the spectral break (as shown in 
Figure [5| as a function of the spectral index obtained when fitting the blazars listed in Table [T] with 



a single power-law. This is plotted in Figure 13(a) A pileup effect would be most pronounced for 



soft sources due to fewer statistics at the highest energies. Consequently, if such an instrumental 
pileup were present, one would expect to see an increase in the measured spectral break as the 
spectral index increases (softens). The best linear fit to the data shown in |13(a) is given by 



Ar(r) = (1. 28 ±0.68)r- (3. 15 ±2. 15), with a fit probability of 22%. This suggests that any pileup 
as a function of spectral index is approximately a 2a effect. 

The next important thing to establish is whether or not this effect could be responsible for 
the ~3o" trend of increasing spectral break with redshift seen in Figure [5j This can be determined 
by looking at the single power-law spectral index distribution as a function of redshift (Figure 



13(b)). In order for a systematic pileup effect to be responsible for the spectral break versus 



redshift distribution observed, the spectral indices of blazars would have to be increasingly soft as 



the source redshift increases. As can be seen in Figure 13(b), this is not the case. The best linear 
fit to the spectral index versus redshift distribution is given by r(z) = (0.46 ±0.87)^; + (3.05 it 0.11) 
and is clearly consistent with being flat. This may seem surprising given that EBL absorption is 
expected to result in softer spectra as source redshifts increase. However, hard spectrum blazars 
are more readily detected by lACTs given that they generally have higher fluxes than their soft 
spectrum counterparts at very high energies. This results in the selection effect that blazars of 
higher redshift are more likely to be intrinsically hard. 

One can conclude from these arguments that a pileup effect does not systematically produce 
a redshift dependent spectral break in lACT spectra. Furthermore, the spectra used in Method 
2 come from three different lACTs all of which use different energy and spectral reconstruction 
techniques. This further reduces the likelihood that a redshift dependent spectral break could be 
purely an instrumental effect. 



- 28 - 



CD 



J3 



il CD 

S ft 

> 

QJ CD 



> 

CD 1i 

P ° 

CO " 



o 

CD 
CD 

o 
o 

CD 



CD 



o 



o 

CD 



t3 

CD 

:3 

03 
CD 

B 

CD 

+^ 

t3 
fl 

03 

> 
CO 



CO 



CD 



CD 
ft 
CD 



.2 

"-+^ 

o 

CD 



" — - CD 
03 3 



i4 CD 



o 

CD 



03 



CD 



o3 
CD 

bC 
03 



o3 



s > 

Co ^ 

iri H 



03 



bO 
O 
iH 

(d 
+J 
ni 
u 



■P 



M 
M 

U 
U 

(d 
\ 

XI 

u 

W 
> 

o 
bO 

nj 
M 
(d 
Pi 

u 

MH 
W 
bO 

•H 

a 
m 



Oh 

O 



03 
O 

03 



O 



CD 



i § 
ft ^ 



03 

> 

CD 

O 



O 
CD 
ft 

> 



CD 

03 



cd 



O 

l-H 



ft 



l-H l-H CO O 

1— I O 00 CO 



00 



lO CM lO ^ Oi 
l-H lO CO CO 00 



>— lOCNi-HOi— "OOi-HOO 

-H-H-H^-H-H-H-H-H-H-H-H 

^ lO CO CD O 

CO CO 02 cn CO 
o o ^ ^ o o o 

I I 



_ . CJ3 CO 00 

^ CO CO CD CM 



COCMCMCOO^COCM-^lOCM^CO 

^^coco^r-cococoi>-m^ 



2 lo S lo ' 

_ ( — ^ --^ . — ^ ' 



o 


o 


o 
O 


o 


CM , 


'St 


CM k 












1 ^ 






"3 










0) 








"fH 


a 


ian 


a 




o 


a 


'a 

o 


"Sb 




o 

Sh 




cS 
















CM CM CM CM CM CM 





00 


CO 




CD 
CO 


CD 
O 


CD 
CM 


CO 




CJ2 
O 


t- 


CO 
CM 


d> 


ci 


d 


d 


d 


d 


d 


d 


d 


d 


d 


d 


-H 


-H 


-H 


■M 


-H 


-H 


-H 


-H 


■M 


-H 


-H 


-H 


lO 


00 
lO 


00 
CM 


o 


lO 

cn 


CM 
CO 


CJl 
CD 


o 

UO 


o 


t- 

CD 


00 
00 


o 


CM 


CM 


oi 


CO 


CM 


CO 




CO 


CM 


CO 





SO 


CD 
O 


CM 

o 


* 

CD 


00 


o 

CM 


o 


* 

CD 
CM 


d 


o 


d 


d 


d 


d 


d 


d 


d 


-H 


-H 




-H 


-H 


-H 


■M 


-H 


-H 




o 


o 


C^i 


o 

CO 




o 

UO 


CD 


i-H 

;D 




CM 














i-H 



00 cri 
CD 
O O 

d d 



^ o 

1^- 00 

o o 









o 


CM 


CD 


00 




CM 


CM 




00 


00 


00 
















d 


d 


d 


d 


d 


d 


d 



5 S 

+ CO 



O CM CJ2 

lO CM 00 

lO CD CO ^ 

+ + 00 lA CM I 

m ^ o !2 

S S o 

CO o - - 



CJ5 

ira 00 

+ 
o 



CM 



CO Cfi 
Dh Oh 



o 

o 



o 
o 

CM CM CO 

1 + + 
+ 00 

CD CM 

O CM CM 

>-3 ^ ~ 



o 



CM 



o 

^ M C« CO 

ffi H H S S 



o 

o "=> 

■ ?j ~r 

» s 

- M " 
S ^ 

^ o 

^ W> • 

S o 

§ S - 

^ ^ 

S o ^ 

^ > g 

■S ^ 

m n t/: 
■« ft 

fe< V fen 



- 29 - 



3- 



UJ 



<l 



-1 - 



-2 



I ' ' ' I ' ' ' I 



1 — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 




LJ I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I J I L. 

0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 

Redshift z 



1. 1ES 2344+514 

2. 1ES 1959+650 

3. PKS 0548-322 

4. PKS 2005-489 

5. RGBJ0152+017 

6. PKS 2155-304 

7. RGB J071 0+591 

8. H 1426+428 

9. 1ES 0229+200 
10.1ES1218+304 
11.1ES 1101-232 
12.1ES 0347-121 



Fig. 5. — Distribution of the observed spectral break versus redshift assuming a break energy of 
1.3 TeV (the choice of which is motivated by the near-IR peak at ~1 //m in the EBL scenarios used 
here). The sample includes the 12 blazars (see Table [T]) available to date with spectra extending 
above and below the break energy. The spectral break is defined as the spectral index below 
the break energy £'b minus the spectral index above the break energy (i.e., AF = F<Eb — r>Eb) 
The black solid line represents the best linear fit given by AF(z) = (8.68 it 5.37)z — (0.24 it 0.71) 
and yields x^/^ = 14.10/10. The grey solid line shows the best linear fit assuming AF(0) = 0, 
giving AT{z) = (6.95 ± 1.65)z with x^/i" = 14.22/11. The long-dashed line shows the AT{z) 
distribution expected from the scenario of Aharonian et al. ( |2006 ), scaled by a factor of 0.55, 
yielding x^/^ = 18.09/12. The short-dashed line shows the AF(z) distribution expected from the 
model of Pranceschini et al. (2008), yielding = 36.23/12. Fitting the data with a flat line at 
AF(z) = yields xV^^ = 31.98/12. 
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Table 2: Summary of the comparison between TeV spectral break data and the EBL models from 
other authors considered here. The columns are as follows: Column 1 : EBL model name. Column 
2: Reference. ColumnS: Best linear fit, of the form AT{z) = mz + b, to the TeV spectral break 
distribution produced by the model (i.e., simulated, not observed data). Column 4- of the 
observed data assuming the best fit from Column 3. Column 5: Confidence level with which the 
model can be excluded given the data. Note: Parenthetical numbers in Columns 4 and 5 are 



obtained using the technique discussed in Section 5.2 



Model 



Reference 



Aharonian ct al. 


(2006 


Aharonian et al. 


(2006 


Aharonian et al. 


(2006 


Stecker et al. 


!006l 


N/A 



Best Fit Ar{z) 



X 



Aharonian 
Aharonian x 0.55 
Aharonian x 0.45 
Stecker (fast) 

Ar(2) = 

Dominguez 
Franceschini 
Kneiske 
Stecker (baseline) 



iDommguez et al.' 
Franceschini et al.. (^2008^ ) 



Kneiske et al. 



Stecker et al. 



(20041 



(p006| 



TAlz - 0.10 
4.05z - 0.05 
3.32^ + 0.00 
-2.832 + 0.42 
0.002 + 0.00 
-0.942 + 0.07 
-1.352 + 0.08 
-I.I82-O.O5 
-4.292 - 0.02 



14.16 (22.95) 
18.09 (21.97) 
19.06 (23.24) 
31.66 (48.54) 

31.98 
34.46 (31.49) 
36.23 (32.67) 
36.46 (34.54) 
61.87 (52.86) 



1.1a (2.2ct) 

1.6(7 (2. la) 

1.7cr (2.2cr) 

3.2cr (4.7cr) 

3.2(7 

3.4(7 (3.2(7) 

3.6(7 (3.3(7) 

3.6(7 (3.5(7) 

5.7(7 (5.1(7) 



10 



20 



30 40 50 



vI^(A-1.6;um) [nW m sr"^] 



Fig. 6. — Constraints, provided by Method 2, on the EBL intensity, in units of nWm~^sr^^, at 
1.6 fj,m and 15 /im. The contours shown are for the 1 (red), 2 (blue), and 3 (green) sigma confidence 
intervals. The intensities for the models of Aharonian et al. ( |2006 ) (red stars), Stecker et al.| ( 2006[ ) 
(green triangles), Franceschini et al. (2008) (blue diamond), Kneiske et al. ( 2004[ ) (orange square), 



and Dommguez et al. (2011) (purple hexagon) are also shown. The three red stars for the scenario 



of Aharonian et al. (2006) represent scaling factors of 1.0, 0.55, and 0.45 (right to left). For the 



model of Stecker et al. (2006), the two green triangles represent the so called fast and baseline 



evolution models (right to left). 
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Fig. 7.— Combined contours from the analysis of RGB J0710+591, lES 1218+304, and Ar{z). 
The contours of lES 1101-232 and lES 0229+200 are left off because they do not further constrain 
the region of allowable parameter space beyond that shown in the figure. The separately colored 
shaded area in each panel indicates the region of overlap for the three contours. 
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Fig. 8.— Combined 2a and 3cr contours from the analysis of RGB J0710+591, lES 1218+304, and 
Ar(z). The contours of lES 1101-232 and lES 0229+200 are left off because they do not further 
constrain the region of allowable parameter space beyond that shown in the figure. Also shown for 
reference are the models of: [Aharonian et ah (2006) (stars) scaled by 1.0, 0.55, and 0.45 (right to 
left); [Stecker et al.| ( |2006'1 ) (t riangles) fast evolution (right) and baseline (left) models; [Franceschini 
et al.[ ( [2008[ ) (diamond) ; [Kneiske et ar[ ( [2004| ) (squ are);| Dommguez et al.[p011[ ) (hexagon). The dot- 
ted line at vl^iVo \im) = 2.0nWm~^ sr~^ indicates the total integrated EBL intensity of Hopwood 
et al. (2010). The dashed lines at vlyj(\h \xyq) = 1.6nWm~^ sr~^ and 1^1^(15 fim) 
indicate the la error on the value. 
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Fig. 9. — Approximate EBL SED 2ct confidence region (blue shaded area) derived from Figure 
[sj Also shown are the models of: Aharonian et al. (2006) (dotted lines) scaled by 1.0, 0.55, and 
0.45; Stecker et al. ([2006 ) (dot-dashed lines) fast evolution (upper) and baseline (lower) models; 
Franceschini et al. (2008j) (dashed line); Kneiske et al. ( 2004| ) (upper solid line); Dommguez et al. 



(2011) (lower solid line and grey shaded area). The filled symbols indicate the galaxy count lower 
limits of (left to right): Gardner et al. (2000) (downward triangles); Madau & Pozzetti (2000) 
(upward triangles); Keenan et al. (2010) (squares with error bars); [Levenson fc Wright] ( |2008| > 



(sideways triangle); Fazio et al. (2004) (squares without error bars); Elbaz et al. (2002) (hexagon 



Hopwood et al. (2010) (circle); Papovich et al. (2004) (diamond). The x at 3.5 /xm indicates the 



tentative detection of Dwek & Arendt (1998). Note: The partially transparent blue shaded area in 
the mid-IR shows the expansion of the 2a confidence region when using the alternate test spectrum 
described in Section 15.21 
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Fig. 10. — Constraints on the EBL intensity, provided by Method 2, using a test spectrum with 6 
spectral points ranging from 400 GeV to 3 TeV. See Figure [6] for a full description of the contours 
and data points shown. 
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Fig. 11. — Figures showing the slope (left) and intercept (right) fit results for the TeV spectral break 
versus redshift distribution (Equation [s]) assuming a particular EBL model and using different size 



error bars for the test blazar spectrum (see Section 3.2). The six different choices of errors shown 



are 15%, 25%, 35%, 45%, 55%, and 65%. The standard deviations (unweighted) of the slope and 
intercept best fit parameters are 0.065 and 0.020, respectively. This spread in values is <3% of the 
error in the best fit determination from observations (Figure^. 
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Fig. 12. — Ar(z) distribution, using EBL absorbed test blazar spectra for a particular EBL scenario, 
calculated after 1, 25, 50, and 100 iterations. 



-37- 




Fig. 13. — (a) Spectral break (as shown in Figure [s]) plotted versus the spectral index obtained 
when fitting the blazar spectra with a single power-law. The best linear fit to the data is given 
by Ar(r) = (1.28 lb 0.68)r — (3.15 ± 2.15). (b) Single power-law spectral index plotted versus the 
blazar redshift. The best linear fit to the data is given by T{z) = (0.46 ± 0.87)^ + (3.05 ± 0.11). 



